Initialisation

## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin  Phi[2] = valeur du noeud  Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
#=======================================#
param <- list(sigma2 = 0.05,
              rho2 = 0.1,
              mu = c(5,90,5),
              omega2 = c(0.5,0.1,0.01) )

#=======================================#
t <- seq(60,120, length.out = 10) #value of times

# --- longitudinal data --- #
set.seed(18031997)
data <- NLME_data(G = 20, ng = 8, time = t, fct = m, param = param)

getDim(data)
##    G   ng    N    n   F. 
##   20    8  160 1600    3
Y <- data$obs

SAEM avec simulation par MCMC

\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]

Phi <- fct_vector(function(sigma2, rho2, mu, omega2) 
  c(-n/(2*sigma2), -N/(2*rho2), -G/(2*omega2), mu), dim = c(1,1,F.,F.) )$eval

S <- fct_vector(function(eta, phi) mean((Y - get_obs(data, eta = eta, phi = phi) )^2 ),
                function(eta, ...) mean(eta^2),
                function(phi, ...) apply(phi^2, 2, mean),
                function(phi, ...) apply(phi  , 2, mean), dim = c(1,1,F.,F.) )

Metropolis Hastings

loglik.phi <- function(x, eta, Phi)
{
  id <- c(1,3,4) #indice de S_1 et S_{3,.} puis S_{4,.}
  sum(Phi%a%1*S$eval(eta = eta, phi = setoffset(x, Phi%a%4), i = 1) + Phi%a%3 * S$eval(phi = x, i = 3) )
}
loglik.eta <- function(x, phi, Phi) 
{ 
  id <- c(1,2)
  sum( Phi%a%id * S$eval(eta = x, phi = phi, i = id) )
}

SAEM

Initialisation

# ---  Initialisation des paramètres --- #
para <- param %>% lapply(function(x) x* runif(1, 1,1.2))
para$rho2 = 0.2 ; para$omega2 <- rep(.1,3)



# --- Initialisation des chaines MC : Z_0 ---
Z <- list(eta = list( matrix(rep(0, G*ng), ncol = 1) ),
          phi = list( matrix(rep(para$mu, G), nrow = F.) %>% t ) )

Etape simulation et maximisation du SEAM

sim <- function(niter, h, Phih, eta, phi, verbatim = F)
{
  M <- length(phi)
  eta <- 1:M %>% lapply( function(i)
    MH_High_Dim_future(niter, eta[[i]], sd.eta , loglik.eta, phi[[i]], Phih, 
                            cores = 1, verbatim = verbatim))
  
  phi <- 1:M %>% lapply( function(i)
    MH_Gibbs_Sampler_future(niter, setoffset(phi[[i]],-Phih%a%4), sd.phi, loglik.phi, eta[[i]], Phih, 
                            cores = 1, verbatim = verbatim )) %>%
    lapply(setoffset, Phih%a%4)
  
  list(eta = eta, phi = phi)
}

maxi <- function(S)
{
  list(sigma2 = S%a%1,
       rho2 =   S%a%2,
       mu =     S%a%4,
       omega2 = S%a%3 - (S%a%4)^2 )
}
sigma2 rho2 mu1 mu2 mu3 omega21 omega22 omega23
Oracle 0.0502965 0.1040354 4.841866 89.83036 5.020998 0.2875514 0.1163807 0.0098201
Initialisation 0.0545425 0.2000000 5.763583 103.74449 5.763583 0.1000000 0.1000000 0.1000000 
niter <- 200
MH.iter <- 10

sd.eta <- .3
sd.phi <- c(.05, .3, .05)

gg <- simulation_test(sim, Phi, param, 100, Z) %>% lapply(function(z)plot(z[[1]], nrow = 2))

res <- SAEM(niter, MH.iter, para, Phi, S$eval, Z, sim, maxi,
            burnin = 100, coef.burnin = 2/3, verbatim = 2)

saveRDS(res, 'saem.rds')
[1] “SAEM execution time = 02min 58sec”
Result of the SAEM-MCMC
sigma2 rho2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3
Real value 0.0503 0.1040 4.8419 89.8304 5.0210 0.2876 0.1164 0.0098
Estimated value 0.0490 0.1016 4.8536 89.8165 5.0815 0.2885 0.1621 0.0345
Rrmse 0.0265 0.0235 0.0024 0.0002 0.0120 0.0034 0.3933 2.5125 
## [1] "SAEM execution time = 02min 58sec"

## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation
## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]